GISS ModelE Description

This document is a short description of what GISS ModelE does and gives some references and descriptions of how it does it. Hopefully this will eventually morph into a full technical paper given enough time and resources!

Model development is an ongoing task. As new physics is introduced, old bugs found and new applications developed, the code is almost continually undergoing minor, and sometimes major, reworking. Thus any fixed description of the model is liable to be out of date the day it is printed. Understanding this, we will endeavour to maintain the web version of the document as closely as we can to the current release, however in is inevitable that some discussion here will on occasion fall behind development.

The GISS GCM model development process is over 25 years old (for a very readable description of the historical development see Hansen et al (2000)). Inevitably, decisions that were made and constraints that existed early on in the process have had influences that are still apparent. While much of the subsequent reworking of the model has led to a reduction in these historical influences, some parts of the model still hark back to the days of punch cards, fortran 66 and line printer output. A charitable interpretation would be that while embracing the new (Fortran 90/95, multi-processing, netcdf, etc.), we endeavour to maintain some of the more harmless GISS traditions (which some might call eccentricities) in a spirit of continuity with those who have previously worked on the model. On the other hand, some of those early decisions (for instance regarding diagnostics, or conservation properties) turned out to be very far-sighted and are a principle reason why the GISS series of models continue to play a useful and important role in the world of GCM simulations. We hope that by continuing to make the GISS models a more accessible and better documented, we will be able to carry on in that vein for at least another 25 years.

CREDITS: This document was mostly written by Gavin Schmidt, with input from Reto Ruedy, Gary Russell, Ye Cheng (PBL physics), M. S. Yao (cloud physics) and Andy Lacis (Radiation). Thanks as always to Jim Hansen and David Rind for supporting the whole ModelE project.

Table of contents

  1. Model files and directory structure

    1. aux

    2. cmor

    3. config

    4. decks

    5. doc

    6. exec

    7. init_cond

    8. model

    9. templates

    10. tests

  2. Overall model structure and main time stepping loop

    1. Initialisation

    2. Main loop

    3. Diagnostics

    4. Input/Output

  3. Atmospheric model

    1. Dynamics

    2. Cloud processes

      1. Moist convection

      2. Large scale condensation

    3. Radiation

    4. Surface fluxes (incl. planetary boundary layer physics)

    5. Turbulence and Dry convection

    6. Stratospheric processes (incl. gravity wave drag)

  4. Sea ice model

    1. Basic thermodynamics

    2. Ice advection

  5. Land Surface model

    1. Soil physics

    2. Snow model

    3. Lake model

    4. Rivers

  6. Ocean models

    1. Imposed Sea surface conditions

    2. Q-flux (mixed layer model)

    3. Russell ocean model

    4. HYCOM (Bleck and Sun)

  7. Tracers

    1. Air mass Tracers

    2. Soluble and Water mass Tracers

    3. Aerosol and Chemistry Tracers

    4. Ocean Tracers

    5. Carbon Cycle Tracers

    6. Special Tracers

  8. References

1. Model files and directory structure

The distribution of modelE has a top directory which contains the following sub-directories: aux, cmor, config, decks, doc, exec, init_cond, model, templates, tests.

modelE/model

This directory contains all the model source code used in all possible configurations. The names of the files are relatively self-explanatory, and each contains groups of subroutines that pertain to a particular function. For instance, SEAICE.f contains modules and subroutines that control the column model for sea ice. The associated file SEAICE_DRV.f contains the interface routines that translate GCM variables to the column model and store the resulting fluxes and diagnostics. Some shared components are in separate sub_directories.

Since the model supports a number of 'plug-and-play' options, there may be a number of files that contain subroutines or entry points (for dummy routines) with the same name. To find out which is relevant for your particular configuration, we recommend automatically generating the documentation for your rundeck using

gmake htmldoc RUN=E001xyz

This will generate a file index.html. From here you can access lists of:

relevant for your rundeck. By descending in to the (highly cross-referenced) tree, you can go from the source file to the routine, to the module whose variables are used to the variable definitions (unfortunately not 100% complete), sideways and back again. (The database variables are the relevant options that you could use in the rundeck to control more aspects of the model's behaviour. NAMELIST variables control the degree of runtime debugging information and diagnostic output and the start and stop dates for the model. These are all fully documented. Trying to work anything out by going directly to the list of all variables is not recommended for the novice!).

To help you get an idea of what is where, here is a brief overview of the source code files:

You will note a couple of modelE conventions. Firstly, declarations of variables and arrays are in modules called *_COM. Column models and local routines do not have a suffix, but the driver routines that mediate between the GCM code and the local physics are call *_DRV.f. Sometimes the _COM declarations are included in the local physics routine. The OPTIONS document explains how the current set of 'plug-and-play' options can be used.

Also in this directory are some template rundecks E001*.R which can be used to generate your own particular configurations. The .depend_E001* are generated by gmake to work out dependencies between source code files so that the proper compilation can be done if one part of the code or rundeck specifications change. Also generated by gmake are the *.sig which are a hack to ensure the time stamps for modules (usually *.mod, but can depend on the compiler) do not lead to circular dependencies.

modelE/decks

This directory should contain all of the user's rundecks, links to run directories, and some useful auxiliary binaries. Since all of these are rundeck dependent, this directory is empty in the distribution (except for the Makefile).

Once you have generated your own rundecks and compiled the model, and/or some auxiliary files, the directory may contain

modelE/exec

The exec directory is where all of the perl and shell scripts that are used are kept. Some of these are called directly from the Makefile and therefore not used directly by a user, but others are used for running, stopping and checking on the run.

modelE/aux

This directory contains a number of auxiliary programs that are useful for the GCM. Many will not be of direct use to most users, but they also give a flavour for how to read and write standard model output.

modelE/doc

As the name suggests, all of the HTML and other documentation can be found here. This includes this files that come with the distribution, but also those that constructed automatically (see above).

modelE/prtdag

Prior to the current incarnation of pdE and POUT_netcdf, it was necessary to have separate post-processing programs to get netcdf format output for the diagnostics. While this is now no longer really necessary, these programs are left in the distribution as guides to how you could produce specialty diagnostics without altering the model code. We do not recommend this in general, but it may be applicable in certain circumstances. Note that keeping these programs up to date has not been a priority, and thus they are not guaranteed to even compile.

2. Overall model structure and main time stepping loop

Initialisation

Initialisation is handled almost exclusively by the call to subroutine INPUT. This routine will initialise all prognostic variables depending on the value of ISTART. Submodules (for instance, ocean, ice, radiation etc.) are initialised through calls to various init_XYZ routines. Most input files will be opened and read at this point.

Main loop

The main loop is contained wholly within the MAIN routine. The order of calling to the physics routines is quite tightly controlled and cannot be changed arbitrarily. Some calls are dummy calls depending on how the model is coupled and which 'plug-and-play' options are in effect. Currently the calling order is: DYNAM (atmospheric dynamics), MELTSI (sea ice lateral melt), CONDSE (moist convection and large scale condensation), RADIA (radiative transfer), PRECIP_XYZ application of precipitation to submodules, SURFCE (calculation of surface heat, water and momentum fluxes), DYNSI (calculation of advective tendencies in sea ice), UNDERICE (calculation of ice-ocean heat, mass and salt fluxes), GROUND_XYZ application of surface fluxes to land surface submodules (earth, land ice, lakes), RIVERF (calculation of river routing), OCEANS (application of surface fluxes and rivers to ocean code), FORM_SI (update of ice module to account for any ice formation in lakes or oceans), ADVSI (advection of sea ice), DISSIP (add heat from atmospheric dissipation to atmosphere), FILTER (pressure filter). At the end of the day calls are made to various daily_XYZ routines that calculate changes (such as in climatologies) that only need be done daily.

Diagnostics

Almost all diagnostics in the model are accumulated and only averaged in the post-processing (or for the end of month printout). Budget page diagnostics are accumulated in the relevant (atmospheric or sea ice) routines, along with relevant atmospheric lat-lon diags. Constant pressure lat-height diagnostics are only calculated in DIAGA, called from within DYNAM. Calculation and print out of the diagnostics is done from print_diags.

Input/Output

Apart from boundary condition files that are controlled by individual modules, restart file input and output is controlled by IORSF. This is a driver routine that calls individual io_XYZ routines for each module. This allows different versions of each module to control the I/O of the required prognostic variables. Depending on the value of iaction in the call, different kinds of I/O can be set (read or write, double or single precision, with or without tracers, with diagnostics or not).

3. Atmospheric model

Dynamics

The solution of the momentum equations is done within the DYNAM. The scheme is leap frog in time with an initial 2/3 time step every 8 leap-frog steps to prevent solution splitting. The dynamics are based on the dry physics (no water vapour effects in the pressure gradient calculation) and use potential temperature as the advected variable. Water vapour and tracers are advected outside the dynamics loop (once every source time step). All temperature, water and tracer advection is done using the quadratic upstream scheme to minimise numerical diffusion.

Cloud processes

CONDSE is a driver that sets up the vertical arrays for the column models for moist convection and large scale condensation, and accumulates diagnostics and output for the radiation and other modules.

Moist convection

The moist convection routine is a plume based model (Yao and Del Genio, 1995) that incorporates entraining and non-entraining plumes, downdrafts (which can also entrain environmental air), subsidence (using the quadratic upstream scheme).

Large scale condensation

The main cloud generating routine LSCOND is based on Del Genio et al 1996, with some modifications to improve the simulation of the nucleation of super-cooled precipitation and the estimate of near-surface cloud formation in very shallow pbl conditions.

Radiation

The radiation code is based on the k-correlation methodology using 33 correlated spectral intervals in the solar and long wave bands (Lacis and Oinas, 1991, Oinas et al, 2001)

Surface fluxes (incl. planetary boundary layer physics)

Surface fluxes are calculated separately for each land surface type (open water, ice, earth and land ice). For each type, a different high resolution PBL calculation is done to extrapolate the first layer atmospheric properites to the surface (defined as 10m above the ground). The diffusion parameters in the PBL layers are functions of the local turbulence closure scheme (Cheng et al, 2002) and are iterated to get a robust estimate of the surface properties (temperature, humidity, and tracer concentration) and fluxes. Generally, 2 or 4 surface flux time steps are done every hour, depending on the vertical resolution.

Turbulence and Dry convection

There are two methods of vertical mixing outside of the surface layer: Either, complete mixing is performed if two boxes are statically unstable (using the virtual potential temperature) (DRYCNV), or a turbulence closure is used to estimate the amount of mixing based on on the grid-scale properties and available TKE (ATURB).

Stratospheric drag (incl. Gravity Wave Drag schemes)

Stratopsheric drag is applied in two ways. Firstly, using a simple formulation as a Rayleigh damping near the top of the model (SDRAG). Secondly, and generally only if the stratsophere is being properly resolved, we have a physically-based estimate of gravity-wave drag detemined from the model simulation of moist convection, mountain waves, shear and deformation (in STRATDYN) (Rind et al, 1988). When gravity waves break, there is also some vertical mixing associated with that. Both GWDRAG and SDRAG are applied within the dynamical leap frog tiome steps to improve stability.

4. Sea ice model

The sea ice model is configured so that it works suitably for the fixed sea ice concentration cases, and also for fully prognostic ocean code. A full calculation is always done for lake ice. Ice advection can be optionally turned on for both dynamic oceans and the Qflux model.

Basic thermodynamics

The sea ice model can be thought of as consisting of two mass layers made up of an upper layer which has a constant thickness (10cm) ice layer and an arbitrary snow amount. The second layer mass has a minimum thickness of 10 cm, but can grow arbitrarily thick. There are two thermal layers in each mass layer which consist of a fixed ratio of the mass. Salinity and tracer values are also calculated on the same mass ratio, with the exception that salinity is assumed to be zero in the snow. At the surface, heat and mass fluxes calculated in the SURFCE routine and the precipitation are applied first to the snow layer and then to the ice. There is advection of heat, salt and mass as a function of the regridding of the mass layers in order to preserve the fixed ratios, i.e. there is a sigma-coordinate for sea ice properties.

The basal fluxes are calculated in the ocean only for the configurations where ice thickness is prognostic (i.e. not for the fixed SST cases). This is done in routine UNDERICE. The fluxes are then applied to both the ice and ocean modules. Similarly, lateral fluxes are calculated in MELT_SI. This is done prior to any surface fluxes calculation to avoid the problem of changing the surface type fractions in the midst of the calculation.

Frazil ice formation can occur in the ocean and lakes in open water conditions and under ice. This ice mass formation rate (and it's associated salt (for the ocean) and energy flux are calculated in either the lake or ocean component and added to the sea ice by routine FORM_SI. There is a minimum lead fraction which is a function of the ice thickness.

Snow can be compacted into ice by accumulating over a maximum thickness, or by rain or surface melting. Surface melting is kept track of since it is used to estimate melt pond formation for the albedo calculation. No specific snow-ice formation is yet included.

Ice advection

Ice advection is a function of the atmosphere-ice and ice-ocean momentum stress and the sea surface gradient. The internal ice pressures and rheology are estimated using the Flato-Hibler viscous-plastic rheology as coded by Zhang and Rothrock, 1999. Due to the requirement that the surface type fractions don't vary over the surface flux calculation we calculate the ice velocities and ice-ocean momentum stress at the beginning of the calculation (in DYNSI), but the ice mass fluxes are not applied until afterwards (in ADVSI). The advection is done using the linear upstream scheme in the two horizontal directions. This code is turned on by using the ICEDYN_DRV ICEDYN modules instead of the ICEDYN_DUM module, and works with any ocean configuration.

In the fixed SST case, using the ice advection does not change the ice variables, but does allow the ice and energy convergences due to advection to be calculated, nad subsequently used in the the q-flux model. The q-fluxes need to be consistent with the ice advection, and so if ICEDYN_DRV ICEDYN is used in the spin-up, it must also be used in the q-flux model (and vice versa).

5. Land Surface model

Soil physics

The soils are modelled separately for bare and vegetated soil. Additionally, there is a canopy layer over the vegetated portion. The soils are modelled using 6 variable depth layers. Evapotransiration takes water from the surface pools, and from below as a function of rooting depths. Runoff is calculated using a TOPmodel approach.

Snow model

The snow model (over land) is a 3-layer snow model (Stieglitz, 1994), that can have a variable extent.

Lake model

The lakes are modelled using a two layer code. The first layer has a minimum depth of 1m, while the second layer is unlimited. Lake mass is effected by precipitation, evporation and lake ice melt and freeze. In addition, river inflow from upstream boxes and outflow to the downstream box modify the first layer properties. There is a small amount of heat diffusion at the interface, and convective overturning occurs if the temperature stratification is unstable.

Rivers

Rivers exist in all land grid boxes and collect runoff from the soil and land ice models. Where there is a lake fraction, the water is considered a lake (as above). Where there is not, there is no intereaction with the atmosphere. Tranpsort of river water depends on the (fixed) downstream direction file and mean slope, and the depth of the upstream river/lake above a predefined sill depth. Some boxes do not contain an outlet, and so can accumulate water indefinitely. As the river water reaches an ocean box, the water is deposited uniformly within that box.

6. Ocean models

Imposed Sea surface conditions

The simplest ocean model is one where the sea surface temperature is read in from a file, either a fixed (annually repeating) climatology, or a transient (monthly varying) realisation based on observations or previous simulations. Monthly means and end of month values are read in and a quadratic approximation used to interpolate to the (fixed) daily value at every grid point. This is constructed to ensure that the monthly mean of the daily interpolated values is equal to the specified mean. Along with SST, the sea ice concentration is also a required field. This is interpolated similarly but is constructed so that the minimum and maximum values are respected. Sea ice thickness is also fixed in these cases, but given the lack of reliable global observations, we set this based on a local scaling to the sea ice concentration.

Q-flux (mixed layer model)

We can calculate the total freshwater mass and heat fluxes into the ice/ocean component from a spin up run with fixed SST. Given a climatology of ocean mixed layer depths, we can calculate the implied ocean heat convergence at the base of the mixed layer. This will incorporate the actual ocean heat transports, but also a residual component related to any errors in the surface heat or mass fluxes. This ocean heat convergence can be input into a thermodynamic mixed layer to give a Q-flux model. Depending on whether ice advection is turned on (see above), the horizontal transport of the sea ice is or is not incorporated into the ocean heat convergence. A full sea ice thermodynamic calculation is performed for this model. Generally this model takes 20 to 30 years to come into thermal equilibrium with any change in forcing, but this can be changed by specifying the maximum mixed layer depth (smaller implies a faster equilibration).

There is an option to allow diffusion into the deep ocean by replacing OCNML with ODEEP in the run deck. This requires a 10 year spin up of the Q-flux model in order to properly set the deeper layers (12 layers, down to 5000m).

GISS Dynamic ocean model (Russell)

This model is a fully dynamic, non-Boussinesq, mass-conserving free surface ocean model (Russell et al, 1995, 2000, Liu et al 2002, 2003). The dynamics is based on a modified Arakawa scheme on the C-grid, with a linear upstream scheme for advecting tracers. Vertical mixing uses the KPP scheme of Large et al (1996). Momentum mixing is modelled as a spatically varying laplacian as in Wajsowicz (1993). The effects of mesoscale eddies and isopycnal diffusion are parameterised as in Gent and McWilliams (1996), but using variable coefficients (Visbeck et al, 1997), and coded as in Griffies (1998).

The model contains up to 12 variable depth subgrid scale straits which contect ocean grid boxes, which would not be connected at the resolution used. In particular, the Straits of Gibraltar, Hormuz, and Nares straits are so modelled. All ocean components are fluxed through these straits as a function of the end to end pressure gradients, balanced against a drag proportional to the straits 'width' which serves as a tuning parameter to get reasonable fluxes.

HYCOM

HYCOM is an isopyncal model based on MICOM that has a hybrid mixed layer scheme near the surface. This has been fully coupled to the modelE atmosphere. For more details contact ssun@giss.nasa.gov.

MOM3

MOM3 is the Modular Ocean Model version 3 from GFDL. For more details on the coupled version please contact ntausnev@giss.nasa.gov.

7. Tracers

ModelE comes with a full complement of tracer code that follows all mass transports in the model. All tracer fields are computed in mass units (not concentration) and thus are perfectly conserved (with one exception explained below). There is a library of pre-coded physics for some common tracers that have been used historically with the GISS models. For the simplest applications, the relevant tracer code needs to be turned on using pre-processor directives in the rundeck (details below), and then by editing TRACER_COM.f the specific tracer can be coded for. Some generic code is allowed for any tracers (such as radioactive decay, diagnostics etc.). Some specific tracers require special physics and these require some extra code (for instance water isotopes or chemistry).

The tracers break naturally into four groups: atmospheric tracers that follow air mass, tracers that are either purely water based or have a soluble component, particulate tracers that may interact with the hydrologic cycle, and tracers that are purely oceanic. These are dealt with in turn below.

Air mass Tracers

The code for tracers of mass transport in the atmosphere is initialised by setting #define TRACERS_ON in the rundeck. These tracers are mostly insoluble gases or ideal tracers such as air mass age or source. All changes in air mass (by advection, convection, the sea level pressure filter) and parameterised mixing (stratospheric diffusion, atmospheric turbulence) are followed by the tracers. There is one exception for the pure 'Air' tracer: this must follow the mass field identically, but during the operation of the sea level pressure filter, a conservation of mass implies a change in concentration. For most tracers, conservation of mass is more important, but for 'Air' concentration is conserved.

Soluble and Water mass Tracers

If the tracer interacts with the hydrologic cycle the preprocessor directive #define TRACERS_WATER should be set in the rundeck in addition to #define TRACERS_ON. This allows for a separate tracer budget in all water reservoirs (cloud liquid water, precipitation, sea ice, lakes, ground water, glaciers and oceans (if a non-trivial ocean model is used)). There can be complicated interactions in the clouds between the in air and in water tracers, and any tracer amount in precip is followed through it's subsequent travels. These tracers are useful for soluble gases and water tracers themselves (isotopes, age, source).

Aerosol Tracers

Ocean Tracers (Russell Ocean only)

These are turned on using #define TRACERS_OCEAN. If no other directives are set, then the tracer code will only be used in the ocean, but if the atmospheric water tracers are also turned on, then tracers can be used throughout the whole system.

Special Tracers

Some special tracer physics come as standard in the model. These relate to tracers that are of special interest to some of the developers.

8. References